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An extended structure-based model based 
on a stochastic eddy-axis evolution equation 

By S. C. Kassinos AND W. C. Reynolds 
1. Motivation and objectives 

Engineering analysis of complex turbulent flows relies heavily on turbulence mod- 
els. A good model should have a viscoelastic character, predicting turbulent stresses 
proportional to the mean strain rate for slow deformations and stresses determined 
by the amount of strain for rapid distortions. Current turbulence models work well 
only in near-equilibrium situations where the turbulent stresses can be predicted 
adequately using eddy viscosity representations. They do not perform well when 
the turbulence is subjected to strong or rapid deformations, which is the case in 
many engineering systems. More elaborate schemes in which the Reynolds Stress 
Transport (RST) equations are included in the PDE system have been used in an 
effort to rectify these problems. While RST models have enjoyed some success, they 
are not yet widely used in industry because they have not proven reliably better 
than simpler models in dealing with the more challenging types of complex flows. 

We have shown that the Reynolds stresses do not always provide a complete 
description of the turbulence state and that this poses a fundamental problem for 
standard RST models that use the Reynolds stress tensor (along perhaps with the 
mean velocity gradient) as the unique tensorial base for the modeling of the unknown 
terms. The inadequacy of component ality information is more pronounced in flows 
with strong mean rotation. These ideas are described in detail by Kassinos and 
Reynolds (1994), hereafter denoted by KR. 

Proper characterization of the state of the turbulence in non-equilibrium flows 
requires the inclusion of structure information to complement the componentality 
carried by the turbulent stresses. We have introduced a number of one-point tur- 
bulent tensors carrying non-local information about the turbulence structure and 
demonstrated how they could be used for the construction of one-point models. 
However, this approach would require the addition of one second-rank and one fully 
symmetric third-rank tensor in the PDE system, a considerable overhead for an 
engineering model. 

These considerations motivated the structure-based model which incorporates the 
key structure information in a simple phenomenological approach. The goal is to 
construct an engineering model with proper viscoelastic character that will reduce 
to the form of a k-e model when the mean deformation is weak, and will match 
rapid distortion theory (RDT) when the mean deformation is strong. 

The backbone of the structure-based model is a one-point, structure-based model 
of RDT for homogeneous turbulence. The development of this RDT model has 
been completed successfully and reported in great in detail in KR, and for that 
reason is not discussed here. This preliminary report focuses on the extensions of 
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the structure-based to flows with weak mean deformation rates. If the structure 
of the turbulence is assumed to be in equilibrium with the mean field and weakly 
anisotropic, the structure-based model reduces to the form of a k-e model. Hence, 
we should be able to extend the model so that it spans between an eddy viscosity 
model, appropriate for weak mean strain rates, and RDT appropriate for high mean 
strain rates. 

2. Accomplishments 

2.1 Overview of the structure -based model 

2.1.1 Algebraic equations 

In a standard k-e model, the turbulent stress tensor Rij is related to the mean 
strain rate tensor S t j through an eddy viscosity 

Rij = p 2 6ij - 2 UrSij Ur = C„k 2 /e (1) 

where Ru = q 2 = 2k. Transport equations are used for k and e but not for Rij itself. 
In the structure-based model, we also determine the Reynolds stresses through 
an algebraic constitutive equation. The difference is that we relate the turbulent 
stresses to parameters of the turbulence structure instead of the mean strain rate: 

Rij = q 2 [(l — <f>)~(Sij — j) + (jxZij + 7 2 -Q^ikyGyj + e jky a yi)\ • (2) 

Here ft* is the mean vorticity vector and fi 2 — The eddy-axis tensor a,j 

carries information on the orientation and shape of large-scale eddies. The two 
scalar parameters (</> and 7) determine the character of the turbulence structure: 
(j> is the fraction of the energy in the jetal mode (motion along the eddy axes), 
1 — (j) is the fraction of the energy in the vortical mode (motion in the plane normal 
to the eddy axes), and 7 is the jet-vortex correlation parameter. In the RDT 
model, we only carry the transport equations for the structure parameters. For 
weak mean deformations (small Sk/e), we need to add the transport equations for 
the turbulence scales k and e. 

The derivation of the algebraic constitutive Eq. (2) for is based on a represen- 
tation the turbulence as a superposition of two-dimensional eddy fields. The moti- 
vation is to account for the effects of the mean deformation on the energy- containing 
eddies. The normalized eddy-axis tensor a\j represents an energy- weighted direction 
cosine tensor of the large eddies, 

A{j = <*«*>) = ^kk^ijf ( 3 ) 

where ( ) denotes an averaging, V 2 is twice the kinetic energy of the basis field, 
and cii is a unit vector aligned with the axis of independence of the field of 2D 
eddies. Note that An = ( V 2 ) = q 2 , and so the eddy-axis tensor scales on the 
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turbulent kinetic energy, as does Rij. When the turbulence structure is isotropic, 
all the eddies are randomly distributed and = Sij/ 3. For turbulence consisting 
of 2D vortices aligned with the x\ direction, an = 1 and all other components 
vanish. Hence, the eddy- axis tensor carries the dimensionality information needed 
by a turbulence model. In fact, a tJ is related to the structure dimensionality tensor 
(see KR) through the model algebraic equation 

Dtj = — a *>)* (4) 


2.1.2 Transport equations 

We use PDE for the structure parameters but not for the turbulent stresses them- 
selves. The evolution equation for the normalized eddy-axis tensor is determined 
from definition (3) and the kinematics of the eddy axis vector a*. In the RDT limit 
a, satisfies the simple equation 


dai 

dt 


— GikO>k G nm a m^n a i 


(5) 


where Gij = Uij is mean gradient tensor. Using (5), definition (3) and some 
analysis, one can show that in the RDT limit the evolution of a,j is given by 


d(%ij 

dt 


+ G* k aki 
O l 


- [3</> + 1 ]Gk m Zkmij + (%<t> 


( 6 ) 


where G*j = Gij — Gkk^ij/3 and S*j = 5,j — Skkbij/3. Closure of (6) requires 
modeling of the energy weighted fourth-moment 

%ijnm = fY Q'iQ'jG'nQ’m) / Q ( 7 ) 


in terms of the second moments a,j. A fully realizable accurate model for Z has 
been developed (see KR). 

The evolution equations for the two scalar parameters 


d(f) __ dj 

dt dt 


(8) 


are derived from the Navier-Stokes equations with some modeling to account for 
information lost in conditionally averaging over the eddies. The exact form of these 
equations is given in KR (see Eqs. 5.10.4 and 5.10.5 therein). 

2.2 Blending of RDT and k-e modeling for homogeneous turbulence 

The simplest approach in extending the RDT structure-based model to slow de- 
formations is the addition of terms in the evolution equations for that model 
the restoration of isotropy as a result of turbulence-turbulence interactions. Similar 
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terms in the evolution equations for <j) and 7 will restore the vortical turbulence 
(<f> = 0,7 — 0 ) appropriate for isotropy. 

For slow deformations we must add equations for the turbulence scales. We 
are currently investigating the use of the familiar k and e equations, with minor 
modifications to take advantage of the structure information provided by the eddy- 
axis tensor. An algebraic relationship is used to obtain the turbulence time scale r 
in terms of k and e. 

The modeling of the return to isotropy terms in the eddy-axis tensor Eq. ( 6 ) is 
perhaps the most sensitive step in implementing these extensions since the simple 
kinematic basis of this equation is critical for full realizability in the RDT limit. The 
added return-to-isotropy terms must capture the key physics without disturbing the 
realizability of the model. For this reason, we next discuss in detail a method of 
extending the a tJ equation that guarantees maintaining realizability. 

2.3 A stochastic eddy-axis evolution equation 

In the RDT limit, the eddy axis vector a, evolves according to the simple kine- 
matic Eq. (5). When the mean deformation is weak, this equation must also involve 
return to isotropy terms accounting for the eddy-eddy (or turbulence-turbulence) 
interactions. Guidance on the form of these isotropization terms can be obtained 
by considering a generalization of the eddy-axis kinematic equation that includes 
stochastic forcing terms, in analogy to the Langevin equation (Arnold, 1974). This 
approach offers the advantage that the realizability of the resulting eddy-axis trans- 
port equation is guaranteed (Durbin and Speziale, 1994). We work with the energy- 
scaled eddy-axis vector 

A, = Fa, (9) 

where V = y/VkVk and a{ is the unit eddy-axis vector. The RDT evolution equation 
for A{ is simply [see (5) and (9)]: 


dAj 

dt 


GikAk G nm a n a m A, G nrn v n v m A, 


( 10 ) 


where Vi = Vi/V. Next we consider a stochastic generalization of (10) given by 


dAi — [GikAk G nm a n a m A, G nm v n v m Ai]dt 

+ C\A{dt + C2dV\?i + CzdWpTpAi + C^eipqdWpAq. 


( 11 ) 


The stochastic forcing in ( 11 ) is provided by the Wiener process dW,(t), which has 
increments that are steps of the random walk and provide Guassian white-noise 
forcing (Arnold 1974). The properties of these increments are 

dWl = 0 dWidWj = dtSij dW t = 0 . ( 12 ) 

The second property in ( 12 ) shows that the Wiener process has magnitude d\V = 
0(dt ) 2 / 2 ; therefore, d\V t /dt is not defined as dt — > 0. Hence, in order to evaluate 
dAij/dt = dA t Aj/dt , we first form the product 


d(A t Aj) — ( Ai + dAi)(Aj + dAj) — A,Aj 


( 13 ) 



Structure-based model 


137 


retaining terms to 0(dt), then average over all eddies, and finally divide by dt. 
Note that the coefficients in (11) are not necessarily constants and are assumed to 
have the appropriate functional forms (in terms of deterministic functions like k , 
e, etc.) that give the correct dimensions in each term. The C\ term is introduced 
by analogy to the Langevin equation and, as will be shown, is needed in order to 
ensure realizability. The C 2 term provides isotropic stochastic forcing that tends 
to randomize the orientation of the eddy sixes. The deterministic vector T, acts 
as an organizing vector for the stochastic forcing in the C3 term; for example, Tj 
could represent an organizing effect for the non-linear interactions provided by the 
structure of the larger scales or the mean field. Finally, the C\ term assumes that 
the non-linear turbulence-turbulence interactions can provide an effective random 
rotation acting on an individual eddy axis. Substituting (11) in (13), one obtains 

d(AiAj} = [GikA^Aj + G jkAfcAi — 2G nm (a n u m ~b VnVm^AiAj + 2C\AiAj]dt 

+ CZdWidWj + C 2 Cz(dW t dW q r q A j + dWpTpAidWj) 

+ C 2 C4(dWi6 jqr dW q Ar + CiptdWpAtdWj) (14) 

+ CZdWpTpAidWqTqAj + C 3 C 4 (e jqr A r Ai + e iqi A t A j )dW q dW p Tp 

~b G^ CtptdVVp-Af €j qr dW? q A r . 

Averaging over the ensemble of eddies and simplifying, one obtains 

^ ~ GikAkj + GjkAki — 2G n m<] 2 Z n mij ~ 2 G nm (v n v m AiAj) + 2C\Aij ^ 15 ) 

+ C^&ij + C 2 T 2 Aij + CzC 4 Tp(tj p tAti + eiptAtj) + Cl(q 2 Sij — Aij) . 

Note that we have no control over the sign of the terms involving C|, C\ and C|, 
which must be positive for realizability, but we have a choice over the sign of the 
terms involving C\ and C3C4. Taking the trace of (15), one finds 

^ = 2P-2e = -2 GnmRnm + 2 C x q 2 + 3 Cl + Cfr 2 ? 2 + 2CW • ( 16 ) 

Therefore we must have 

2 Ciq 2 + 3 C\ + C\Y 2 q 2 + 2 c 2 q 2 = -2 c . (17) 

Based on dimensional considerations, we let 

Ci = Ci/t C 2 = C 2 V~e C 3 = C 3 /V7 C 4 = C 4 /^. (18) 

Then condition (17) becomes 


2 Ciq 2 /T + 3 C 2 e + C 3 r 2 q 2 /T + 2 C\q 2 jr = -2e , 


(19) 
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and if we assume r = q 2 / 2e then 


2 Ci + \cl + cl r 2 + 2 cl = -l . 


( 20 ) 


Note that realizability requires C\ < —1/2. Using (3), (15), and (18), one finds that 
the evolution equation for is given by 


da 




dt 


— G{ k&kj ~\~GjkO>ki *^G nrn Z nrn ij ^G nm {v n V m AiAj) j q 2S nm r nm a 




+ 1 [(1 + 2Ci + ClT 2 - Cl) aij + l -{Cl + 2C\)6,j + C 3 C 4 r p (e jpt a ti + e ipt a tj ) ] . 

( 21 ) 

Note that the trace condition dau/dt = 0 is satisfied if (20) is assumed. Next we 
require that the terms involving Ci, C%, Cf, and C\ in (21) combine to produce a 
return to anisotropy term of the form Cl(\6ij — aij). This can be accomplished if 


3 

2< 


-1 = 2 C\ + ^ C I + C 3 2 r 2 + 2 C 2 . 


( 22 ) 


Note that (22) is identical to (20), and hence this additional requirement is auto- 
matically satisfied if the trace condition (20) is imposed. Solving (20) for C\ and 
substituting back into (21), one finds 


da 


dt 


, — (*ik a kj "1“ G jk a ki 2 G nrn Z nrn ij 2,G nm (v n v m AiAj)/q + 2 S nrn V nm aij 


1 

+ - 
r 


2^2 + )(“^tjf — aij) + CsCtTp^tjptati + * 


(23) 

The fourth term on the RHS of (23) can be evaluated using the conditional averaging 
procedure described in KR (see pp. 85-95). Substituting the resulting expression 
in (23) produces 


da 


dt 


, ~ G%k°kj + Gjk a ki [3<£ + \\G k m Zkmij + (3<^ 1 )G nm a nm aij 


41 ^Jt 

27*$n m ~Q~tnkt(Ztmij Q'tm a ij) 


(24) 


H — ^(^2 + ^ C 2 ){-r8ij — aij) + CzCiTpfaptati + eip t a t j) . 
T H O J 


The algebraic k-e Eq. (1), expressed in terms of the Reynolds stress anisotropy 
fij = ri j - Sij/ 3, is given by 


2vt (T 

r — —L q* — —or — q* — ~c _ q* 

lJ - q 2 ij ~ 'Ue ” M 2e 0 


(25) 
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where vt — C^k 2 /e = Cf,(q 2 ) 2 /( 4e) is the turbulent viscosity. For irrotational mean 
strain, the algebraic constitutive Eq. (2) for the structure-based model produces 
(with <f> = 7 = 0) 

1 . 




- 2 a H 


(26) 


where dij = aij — \6{j is the anisotropy of a, j . FVom (25) and (26) we see that 
consistency with k-e modeling (in the weak strain limit) would require 




Now if we assume equilibrium under weak strain rates in (23), we obtain 

l[l(C 2 2 +2C 2 4 )a ij ] = ±S: j . 


(27) 


(28) 


Substituting (27) in (28), one finds that consistency between the structure-based 
model and k-e modeling in the limit of equilibrium structure under weak deformation 
requires 


l(C 2 + 2C 2 )= 4 


15 C u 


(29) 


Next we consider two limiting cases where there are no non-linear eddy-eddy inter- 
actions, and hence the coefficient §(C2 + 2 C\) should vanish. The first case is that 
of a 2C-field of jets having a 2 2 = 0 and 0 = 1, corresponding to the type of structure 
one might expect to find at the wall in a boundary layer. The jets in this 2C field 
have no way of re-orienting each other towards a more isotropic distribution. The 
second limiting case is that of a 2D-field of vortices with an = 1 and 0 = 0, cor- 
responding to the RDT limiting state in the irrotational axisymmetric contraction 
flow. Again these vortices have no means of re-orienting each other, and the return 
to isotropy must shut off. Both of these limiting cases can be accounted for by the 
postulating the functional form 


-(<?!+ 2 cj) = 0(1 - m -- 2 ) 


(30) 


where a 2 = a,* a*,. Then the k-e consistency requirement (29) for equilibrium under 
weak strain rates produces 

2 


a — 


5 C u ‘ 


Based on this analysis, we propose using 


da 


dt 


_ 


— Gik a kj + G* k aki - [30 + l)Gl m Zkmij + (30 - 1 )G£ m a nm aij 

* 

— ^7 Snm^Q e nkt(Ztmij ~ ^tm^ij) 

+ ~ [ a (l ~ ^)(l — a2 )(^ij ~ °ij) + C e -^-(ej pt au + Ciptat,)j • 


( 31 ) 


( 32 ) 
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Considering equilibrium in homogeneous shear suggests the values using a « 1.8 
and C e « —0.35, which we have adopted. Note that we have taken the deterministic 
vector r t to be a unit vector aligned with the mean vorticity vector ft*; this choice 
was suggested by looking at a number of homogeneous flows, including homogeneous 
shear, irrotational axisymmetric strain, and plane stain. 

2.4 Extension of the scalar transport equations 

We are currently using one simple term in each of the two scalar equations that 
tend to restore vortical turbulence (<j> = 0,7 = 0), appropriate for isotropy. The 
form of the extended equations is as follows: 


dt 


RDT ... - C+<f>/r 


with C<f> = 1.3 


(33) 


d'y 

dt 


RDT ... - C 77 /t 


with C 1 = 2.8 . 


(34) 


The numerical values for the model constants were calibrated for homogeneous 
shear. Here RDT stands for the RHS of these equations in the RDT limit as given 
in KR [see Eqs. (5.10.4) and (5.10.5)]. 


2.5 Evolution of the turbulence scales 

The choice of turbulence scales to be used in a turbulence model is not unique. For 
example, standard k-e models use transport equations for k and e, and determine 
the turbulence time scale through an algebraic equation. Another possibility is to 
use the evolution equation for the time scale r (or the reciprocal time scale w) along 
with the equation for fc, and then evaluate e from an algebraic equation. Each of 
these approaches has some problems. We are currently investigating the use of the 
k and e equations in the form shown below. 


t= p -‘ < 35 > 

% = l-c.s» - C.Syry - C 0 /r - C a (36) 

at 

with 

C d = \ C s = 3.0 Co = — Cu= 0.01 . (37) 

3 o 

Note that the e equation has the standard form except for the last term involving 
Cn • This term takes advantage of the structure information in a x] and allows for a 
decrease in the dissipation rate in the presence of mean rotation, except when the 
turbulence becomes two-dimensional, as observed in direct numerical simulations. 


3. Evaluation of the proposed extensions 

In this section, the extended structure-model given by (2), (4), and (32)-(37) 
is tested for four independent homogeneous flows. First we summarize the values 
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of the constants that we will be using. For clarity we also include the values for 
the constants in the RDT model [Eqs. (5.9.6)-(5.9.11) in KR] for which we use the 
notation of KR: 


• Rapid model : 

Ct = 5.9 

C 2 = 2.0 

C 3 = 7.0 C 4 = 2.5 

(38) 

• Slow model : 

a = 1.8 

C'e = -0.35 

C+ = 1.3 C 1 = 2.8 



*-i 

C s = 3.0 

Co = — Cn = 0.01 . 

0 

(39) 


3.1 Homogeneous shear in a rotating frame 
We first consider the problem of homogeneous shear in a rotating frame. The 
mean velocity gradient tensor Gij, the frame vorticity £l{ , and frame rotation rate 
Cti are defined by 

/0 5 0\ 

= I o o o , 2 ft, = ft/ = (o,o,n f ). < 40 ) 

\0 0 0 / 

We consider initially isotropic turbulence 

r ij = = Jtoi — € 0 . (41) 

First, we consider the case of homogeneous shear in a stationary frame (fif = 0) with 
an initial Sk o/eo — 2.36. Figure 1 shows the model predictions for the components 
of the normalized Reynolds stress tensor = Rij/q 2 . The symbols are from the 
direct numerical simulation of Rogers et al (1986), which also had Sk 0 /e 0 = 2.36. 
The agreement between the model predictions and the direct numerical simulation 
is good. As shown in Table 1, the equilibrium state predicted by the model is in 
good agreement with the experiments of Tavoularis k Karnik (1989). 



FIGURE 1. Time evolution of the normalized Reynolds stress tensor in homoge- 
neous shear for Sko/eo = 2.36. Comparison of the predictions of the structure-based 
model ( ) with the direct numerical simulations of Rogers et al. (1986) (• ). 
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The solution in the case of homogeneous shear in a rotating frame depends on the 
initial conditions only through the dimensionless parameter SfcoAo, and on the 
frame vorticity through the dimensionless parameter Q.? / S (Speziale et. al, 1991). 
The value of & /S determines whether the flow is stable, in which case k and e 
decay in time, or unstable, in which case both k and e grow exponentially in time. 
In the stable regime (e/Sfc)oo = 0, and in the unstable regime (e/Sk ) > 0. 


Equilibrium 

Values 

Structure 

Model 

Experiments 

rn 

0.53 

0.51 ± 0.04 

r 22 

0.18 

0.22 ± 0.02 

T 33 

0.29 

0.27 ± 0.03 

r 12 

-0.16 

-0.16 ± 0.01 

Sk/e 

5.30 

4.60 ± 0.50 

P/e 

1.70 

1.47 ± 0.14 


TABLE 1. Equilibrium results for homogeneous shear: comparison with the exper- 
iments of Tavoularis & Karnik (1989). 



FIGURE 2. Bifurcation diagram of the structure-based model for homogeneous 
shear in a rotating frame. 

Linear analysis and LES show that the flow is unstable for —0.21 ^ & / S < 1 and 
stable outside these bounds. The most unstable case, having the highest growth 
rate for k and e and the largest (e/Sk ) <», corresponds to f l f / S = 0.5. Figure 2 
shows the bifurcation diagram for the structure-based model. The structure-based 
model does an excellent job predicting the location of the bifurcation points A and 
B, and that of the most energetic state C (largest growth rate for k). 

In the absence of DNS or experimental data, we evaluate the model performance 
using the large-eddy simulations of Bardina et. al (1983). Figure 3 shows the 
evolution of the normalized kinetic energy k/k 0 with non-dimensional time St. Note 
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FIGURE 3. Time evolution of the turbulent kinetic energy in rotating shear flows. 
Comparison of the predictions of the structure-based model (lines) with the large- 

eddy simulations of Bardina et. al (1983) shown as symbols: Sl f = 0, ( , ■); 

nf = 0.5 5, ( , • ); (If = S, ( , 4 ); & = -5, ( , ♦). 

that the model captures the general trends correctly. For example, it correctly 
predicts that the highest rate of growth (for both k and e) should occur for SI* = 5/2, 
which RDT shows is the most unstable case. It also predicts a weak rate of growth 
for the case = 5 and a decay (relaminarization) for & = —5. The numerical 
agreement with the LES is reasonable, but the model tends to predict somewhat 
lower rates of growth, particularly so in the case SI? = 0.55. This problem is 
common to all the currently available second-order closures as noted by Speziale et. 
al (1989). 

S.2 Axisymmetric strain 

Next, we consider the performance of the extended structure-based model for the 
cases of axisymmetric contraction and expansion in homogeneous turbulence. The 
mean velocity gradient tensor is given by 


5 0 0 

Sij = | 0 -5/2 0 

0 0 -5/2 , 


(42) 


with 5 > 0 for contraction and 5 < 0 for expansion. We consider an initially 
isotropic state as specified in (41). The solution depends on these conditions through 
the non-dimensional parameter 5fco/c 0 . Comparisons are made with the DNS of 
Lee & Reynolds (1985). In both cases, we compare with the slowest runs from 
these simulations, which correspond to Sko/to = 0.56 (contraction case AXK) and 
Sko/e 0 = 0.41 (expansion case EXO). 

In Fig. 4(a), we consider the time evolution of the components of the Reynolds 
stress anisotropy r^. The total strain 




( 43 ) 
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FIGURE 4. Comparison of the model predictions ( ) with the direct numerical 

simulations of Lee & Reynolds (1985) (• ) for irrotational axisymmetric contraction 
with Sko/eo = 0.56. (a) Evolution of the Reynolds stress anisotropy tensor fij. (b) 
Evolution of the normalized turbulent kinetic energy k/k 0 . 



C C 


FIGURE 5. Comparison of the model predictions ( ) with the direct numerical 

simulations of Lee & Reynolds (1985) (• ) for irrotational axisymmetric expansion 
with Sko/eo = 0.41. (a) Evolution of the Reynolds stress anisotropy tensor f tJ . (b) 
Evolution of the normalized turbulent kinetic energy k/ko . 



C C 

FIGURE 6. Comparison of the model predictions ( ) with the direct numerical 

simulations of Lee & Reynolds (1985) (• ) for irrotational plane strain with Sko/e 0 = 
0.50. (a) Evolution of the Reynolds stress anisotropy tensor f, r (b) Evolution of 
the normalized turbulent kinetic energy k/ko ■ 
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serves as the time coordinate in this plot. The model predictions, shown in solid 
lines, are in very good agreement with the DNS of Lee & Reynolds (1985) shown 
in symbols. The same good agreement between model and DNS is obtained in 
Fig. 4(b), where we consider the time evolution of the normalized turbulent kinetic 
energy k/ko. 

The evolution of the Reynolds stress anisotropy in the case of axisymmetric ex- 
pansion is considered in Fig. 5(a). Note that the model underpredicts significantly 
the level of anisotropy as compared to the DNS of Lee &; Reynolds (1985). The 
model prediction for the evolution of the normalized turbulent kinetic energy [shown 
in Fig. 5(b)] is accurate up to C « 2.6, but eventually it also degrades. This defi- 
ciency of structure-based model is discussed shortly. 

S.3 Plane strain 

We now turn to the case of homogeneous turbulence subjected to plane strain. 
The mean velocity gradient tensor is given by 

(S 0 0\ 

Sij= 0 — S 0 . (44) 

\0 0 0 / 

We consider initially isotropic conditions corresponding to (41) with Sk 0 /e o = 0.5. 
These conditions correspond to the slowest run (case PXA) reported by Lee L 
Reynolds (1985). 

Figure 6 (a) shows the time evolution of the Reynolds stress anisotropy. Note 
that the model predictions are accurate only for very small total strain and quickly 
degrade, particularly for rji and r 33 . As in the axisymmetric expansion case, the 
model prediction for the rate of decay of k/ko remains accurate for a somewhat 
larger total strain, but eventually it degrades also [see Fig. 6 (b)]. 

8.4 Some problems with the current approach 

The relatively poor performance of the structure-based model in the axisymmetric 
expansion and plane strain flows prompted us to take a closer look at both the 
physics of these flows and at our model. What we have learned helped us understand 
better these flows and also provided us with a solution to the problems faced by the 
structure-based model in these flows. 

Rapid distortion analysis (RDT) shows that under irrotational mean strain fij = 
dij. This result is clearly exhibited in the most rapid runs from the DNS of Lee & 
Reynolds (1985), including the rapid expansion and plane strain runs corresponding 
to Sfco/eo = 41.0 and Sko/eo = 50.0 respectively. However, when the slowest 
runs for these two flows are considered, corresponding to cases EXO and PXA 
discussed above, one finds that ? d (see Fig. 7). These observations become 
even more interesting if one notices that the level of stress anisotropy r,j in the slow 
axisymmetric expansion and plane strain runs exceeds the level of stress anisotropy 
in the corresponding rapid runs! This effect is demonstrated in Fig. 8 , where we 
show plots of \II r \ — fijfji/2 versus |JJd| = dijdji/ 2 . The open symbols correspond 
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Figure 7. Comparison of the Reynolds stress anisotropy r t] ( ) with the 

dimensionality anisotropy d,j ( ) from the direct numerical simulations of Lee 

& Reynolds (1985). (a) Irrotational axisymmetric expansion with Sk 0 /e 0 = 0.41. 
(b) Plane strain with Sko/eo = 0.50. 
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FIGURE 8. The second invariants of the stress anisotropy — II r vs. the second 
invariant of the dimensionality anisotropy — lid- (a) Axisymmetric expansion at 
Sko/eo = 0.41 (■) and at Sko/eo = 41.0 (0). (b) Plane strain at Sko/eo = 0.50 (■) 
and at Sko/eo = 50.0 (0). 

to the most rapid run and the closed symbols to the slowest run of Lee &: Reynolds 
(1985) for each flow. Note that in the rapid runs \IId\ « \II r \ whereas in the slow 
runs \IId\ < |IJr|. What is more, in each flow, the maximum level reached by \II r \ 
is higher in the slow run that it is the rapid run. 

By using a linearized version of the RDT evolution equations for f{j and J tj , valid 
for small anisotropies, we have been able to show that these intriguing effects are 
primarily controlled by the rapid terms in the two evolution equations. In other 
words, RDT will maintain r = d if it is initially true, but an arbitrarily small 
deviation A = r — d will be amplified by the rapid terms. The initial conditions of 
the simulations of Lee Sz Reynolds imposed a very small initial Ao = ?(0) — d(0). 
However, even in the absence of any initial A 0 , such a deviation could be triggered 
by unequal rates of return to isotropy for the two tensors. 

The fact that these unexpected effects (once triggered by the initial conditions or 
non-linear effects) seem to be dominated by the rapid terms prompted us to take 
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a closer look at our rapid model. We believe that the difficulties encountered in 
these flows are related not to the slow model developed above, but rather to the 
form of the basic constitutive Eq. (2), which relates the Reynolds stresses to the 
eddy-axis tensor. The reason for this failure lies in the fact the current version of 
the structure model assumes that, in absence of mean rotation, <£ = 7 = 0 . This 
means that the principal axes of r x j remain locked onto the principal axes of a x j. 
This is appropriate for the RDT of initially isotropic turbulence, when the eddies 
do not have time to interact with each other. The cases examined above show that 
this is not appropriate for slower mean strain rates where the non-linear eddy-eddy 
interactions are important. These non-linear eddy-eddy interactions provide an 
effective eddy rotation acting on an individual eddy due to the circulation associated 
with the background sea of eddies. The effective eddy rotation tends to rotate the 
principal axes of the stresses associated with an individual eddy so that these become 
misaligned with the eddy axis, and some <j> and 7 are produced. But in order to 
capture these effects it is not enough to allow for non-zero </> and 7 under irrotational 
strain; we also need to replace the mean vorticity ft; in ( 2 ) with the effective eddy 
rotation rate ft*; this will produce a contribution in the jet- vortex correlation term 
even in the absence of mean rotation. Simple kinematic analysis (see Appendix I 
in KR) shows that ft* is given by 


ftj — tirpQ'pQ'r 


a i = 


dcii 

dt 


(45) 


Note that because of (45) the effective eddy rotation rate will be sensitive to the 
slow model adopted in the a x (or A*) evolution equation. Some preliminary analysis 
suggests that these changes in the constitutive Eq. ( 2 ), coupled with an appropriate 
slow model in the A, (and hence a,j) equations, will allow the structure-based 
model to access states above the RDT limit on the axisymmetric expansion line of 
the anisotropy invariant map. 


4. Summary and future plans 

We have proposed and implemented an extension of the structure-based model 
for weak deformations. It was shown that the extended model will correctly reduce 
to the form of standard k-e models for the case of equilibrium under weak mean 
strain. The realizability of the extended model is guaranteed by the method of its 
construction. The predictions of the proposed model were very good for rotating 
homogeneous shear flows and for irrotational axisymmetric contraction, but were 
seriously deficient in the case of plane strain and axisymmetric expansion. 

We have concluded that the problem behind these difficulties lies in the algebraic 
constitutive equation relating the Reynolds stresses to the structure parameters 
rather than in the slow model developed here. In its present form, this equation 
assumes that under irrotational strain the principal axes of the Reynolds stresses 
remain locked onto those of the eddy-axis tensor. This is correct in the RDT 
limit, but inappropriate under weaker mean strains, when the non-linear eddy-eddy 
interactions tend to misalign the two sets of principal axes and create some non-zero 
<f> and 7 . 
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We plan to modify the constitutive equation and the evolution equation for the 
eddy- axis tensor as necessary to reflect these effects. This will require replac- 

ing the mean vorticity vector ft, in the constitutive equation by an effective eddy 
rotation rate fi* = e t > p a p d r that correctly accounts for the non-linear effects de- 
scribed above. The slow model in the eddy-axis equation may have to be adjusted 
accordingly since the effective eddy rotation rate fi* will be sensitive to it. 

Once these modifications have been implemented and evaluated, we will focus in 
extending the structure-based model for inhomogeneous flows. This extension will 
require the addition of diffusion terms in the transport equations for the structure 
parameters and the turbulence scales. Some preliminary work in determining the 
form of the diffusion terms and appropriate boundary conditions for these equations 
has been carried out. 
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